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Summary 

A computer program based on an integral equation due to Hufford has been 
written to calculate groundwave field strength at low and medium frequencies. Account 
is taken of arbitrary variation of the ground constants and the height above sea level 
along the path of propagation. A variety of tests show the program to be fast and the 
results accurate. 

The results confirm the widely held belief that departure of the surface from a 
smooth sphere is not usually very important in terrain that is not mountainous. Never- 
theless a useful gain in groundwave field strength can result from siting a transmitter on 
high ground. 
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List of principal symbols 

Rationalised units are used and the time factor 
exp(jotf) is suppressed. Where no unit is given below the 
quantity is dimensionless. The units are not applicable to 
the FORTRAN data specified in Section 8.3. 



A,B,P 

a 

B 

CUVJ) 

D 

E(ND,ID) 



e A ,h A ,e B ,h B 



F 
f(m) 

G = Z A b/^A B 



/ 

K(N,I) 

k 

N 

n 

p 1 = pexp(-jb) 

R (=ND) 
r (=ID) 

Z AB 

277 

e r = jfc-/1-8x 10 10 ct/F 
V 



points defined in Fig. 1 
radius of Earth (m) 
constant defined by Equation 
(10) (n\- % ) 
see Equation (9) 
calculation interval (m) 
effective surface impedance de- 
fined by Equation (11) (ohms) 
electric and magnetic fields per 
ampere impressed between the 
terminals of A(B) (ohms/m, 
m" 1 ) 

values of e A ,h A ,e' B ,h' B when the 
ground is plane and perfectly- 
conducting 
frequency (Hz) 

function of integer as defined by 
Equation (13) 

complex groundwave attenua- 
tion factor 

number of intervals into which 
r is divided 

V-i 

coefficient defined by Equation 
(14) 

dielectric constant (i.e. real part 
of relative permittivity) 
number of intervals into which 
R is divided 

unit vector normal to the surface 
Norton's complex numerical dis- 
tance (defined by Equation (19)) 
range defined by Equation (5) 
(m) 

factor multiplying Earth's radius 
to take account of refraction 
distance defined by Equation (5) 
(m) 

mutual impedance between two 
vertical grounded aerials (ohms) 
value of Z AB when the Earth is 
plane and perfectly-conducting 
(ohms) 

free space propagation constant 
(radians/m) 

complex relative permittivity 
surface impedance (ohms) 
intrinsic impedance of free space 
(1207rohms) 



P.Pa-Pb 
a 



free-space wavelength (m) 

phase retardation defined by 

Equation (5) (radians) 

radial distances defined in Fig. 1 

(m) 

conductivity (S/m) 

angles defined in Fig. 1 (radians) 



1. Introduction 

The planning of l.f. and m.f. broadcast services entails 
the computation of field strength over non-uniform ground. 
This is required both for the assessment of service areas and 
for the determination of the radiation pattern and efficiency 
of an aerial, with its Earth system, from measured field 
strengths. 

Groundwave propagation over imperfectly conduct- 
ing Earth has been the subject of theoretical research since 
the beginning of the century and is still receiving attention; 

1 2 

the reader is referred to articles by Wait ' for a review of 
published work and for self-contained derivations of the 
principal results. Even for a homogeneous plane Earth 
Sommerfeld's original solution is not simple; it is equiva- 
lent to resolving spherical waves diverging from the source 
into an angular spectrum of plane waves, and supposing 
each elementary plane wave to be reflected by the ground 
in accordance with the appropriate Fesnel reflection 
coefficient. 

In order to take the spherical form of the Earth into 
account, the continuous angular spectrum of plane waves is 
replaced by a series of zonal harmonics, which is slowly 
convergent but has been transformed into a more useful 
series by Watson. A thorough exposition of the theory 
has been given by Bremmer. 

The results of the classical theories referred to above 
were reduced to a convenient form for engineering use by 
K.A. Norton. 4 ' 5 ' 6 

The effect of inhomogeneous ground has been studied 
by many workers from 1944 onwards, using three principal 
methods. The first of these, due to Feinberg, is a pertur- 
bation method based on Green's theorem.* Secondly, 
p 

Clemmow used dual integral equations to study propaga- 
tion across a straight boundary between regions of a plane 
Earth having different ground constants; this elegant 
method cannot readily be extended to more complex cases. 
The third approach, which is the basis of the program to be 
described, is a perturbation method based on the compen- 
sation theorem, which leads to substantially the same 
mathematical formulation as does the Green's theorem 
method, but seems more easily understood. Using a com- 

* See Reference 7, Chapter 9 for references and discussion. 
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puter it permits both non-uniform surface impedance and 
non-uniform curvature to be taken into account, and pro- 
vides a convenient treatment of a uniform plane or 
spherical Earth. 

The radius of curvature of a smooth spherical Earth 
has a considerable effect on the field strength at longer 
ranges but it has not in the past been customary to take 
into account departures of the Earth's surface from a 
smooth sphere. Although experience has shown that the 
effect of this is usually too small to warrant consideration 
when computing service areas, it was thought possible that 
local irregularities could have an influence on the compu- 
tation of effective radiated power, and of conductivity, 
from measured field strengths. Moreover it could some- 
times be desirable to take account of ground contours in 
choosing the site of a high-power transmitting station. 
When the transmitter generates several hundred kilowatts 
even a 1 dB increase in field strength is worth a considerable 
annual sum. 

it was decided to develop a program aimed at the 
following objectives: 

(i) To provide means of determining the effect of ground 
contours in the vicinity of the transmitter when there 
is particular reason to suspect that this could be 
important. 

(ii) To determine whether the effect of the ground con- 
tours is sufficient to be worth taking into account in 
routine planning. 

(iii) To provide an accurate means of checking cheaper 
(i.e. faster) methods of calculation that might be in- 
corporated in a program for routine service planning. 




centre of Earth 



Fig. 1 - Section of path 



plane and perfectly conducting (the broken line AB). 
The object will be to determine the complex attenuation 
factor G given by 



G-Z AB /Z AB 



(1) 



Suppose that initially the Earth is plane and perfectly- 
conducting and that unit current impressed between 
terminals A results in electric and magnetic fields* e A , h B . 
Similarly, if unit current were impressed between terminals 
B, electric and magnetic fields e B , h B would result at A. 
Now suppose that the Earth becomes curved and imper- 
fectly-conducting and that in consequence e A , h A , e B , h B 
change to e A , h A , e B and h B . At the same time Z AB 
changes to Z AB . Then one formulation of the compen- 
sation theorem is** 



(iv) Although it was expected that the program would be 
too expensive in running time for use in routine plan- 
ning it seemed just possible that it could ultimately 
be refined to make it usable in this way. 



Z' 



AB 



: ^AB + 



xh A -e A xh B 



.ndS 



(2) 



2. The integral equation 

The integral equation is essentially that derived by 
Hufford 10 using the Green's theorem method referred to 
in the Introduction. However here we will follow Wait and 
others* in using the compensation theorem, which is more 
readily understood by many engineers. The analysis is to 
be found in Reference 7, Chapter 8, but will be summarised 
for convenience here. 

In Fig. 1 the solid curve line AB represents the surface 
of the Earth. A and B are two aerials; of these A will be 
regarded as the transmitting aerial, while B is regarded as a 
receiving aerial used to determine the field strengths in 
terms of the mutual impedance between A and B. Let 
Z AB be the (complex) value of this mutual impedance, and 
let Z AB be the value that it would have if the Earth were 



' See Reference 2 for a useful survey and further references. 



where the surface S includes the region of the Earth's 
surface making a significant contribution to the integral. 

Equation (2) is simplified as follows: 

(i) e B , h B and Z AB are evaluated for a perfectly-con- 
ducting plane Earth. 

(ii) e A' ^A and Z'j^q are expressed in terms of the 
attenuation factor G. 

(iii) It is assumed that the relationship between the com- 
ponents of the electric and magnetic fields E and H 
tangential to the Earth's surface may be expressed in 
terms of a (complex) surface impedance 7? as fol- 



7,9 



lows 



* Lower-case symbols denote fields due to unit impressed current. 
** Reference 9, Equation (4), or Reference 7, Equation (8.25). 
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E = i?n x H 



(3) 



where n is a unit vector in the direction of the upward 
normal. 

(iv) The gradient of the surface, relative to the tangent 
plane at A, is small. 

Equation (2) becomes:* 



where e r , the complex relative permittivity, is given by 

e r = k-j l-8x10 lo a/F (8) 

where k is the dielectric constant, 

a is the conductivity in Siemens/m, 

F is the frequency in Hz. 



G(p<J>) = 1 



2?r 



* + ■ 



-//3 [(AP) + (PB)-(ABj] 



where j3 Q = 27r/A, X being the free-space wavelength and rj 
is the intrinsic impedance of free space. 

Next the double integral is reduced to a single integral 
by integrating with respect to a co-ordinate directed along 
the surface perpendicularly to the path of propagation. 
Variations of the surface impedance and the height along 
the perpendicular are assumed to be sufficiently slow so that 
the principle of stationary phase permits them to be ignored. 
It is also permissible to approximate an element of length 
measured along the surface by ad§, and to replace (AB), 
(AP), and (PB) by distances measured along the mean 
sphere, except in the exponent. A further approximation is 
to neglect reflections due to irregularity of the ground and 
to discontinuities in its surface impedance. For example, 
there is negligible reflection back towards the transmitter 
by a land/sea boundary. This assumption, which is sup- 
ported by experience, enables the range of integration to be 
restricted to the line between transmitter and receiver. 
Writing 

| =j8 [(AP) + (PB)-(AB)] 
R=a§ 



(5) 




G(r)dr 



(In the special case of a plane Earth, \p and % are zero.) 
The relative surface impedance 7?/t? is given by* 



(6) 



— ~(e r + 1)- 



(7) 



Reference 7, pp. 114 and 1 15. 
* Reference 7, Equation (6.8). 



(AB) 
(APMPB) 



G(p<p)dS 



(4) 



3. Method of so! tits on 

Equation (6) expresses the attenuation factor at a 
range R in terms of its value at all ranges between zero and 
R. It is solved in a progressive manner; once G(r) has 
been derived for a series of equally-spaced values of r 



r=Q,D,2D, 



,nD 



one more application of the equation gives a new value 
G([n+\]D). 

The first step is to replace the integral by a sum, 
replacing 7? by ND and r by TD to give 



TV 



"5} 



GiND) = 1 -BD /2 J E{ND,ID)C(N,I)G{ID) (9) 

/= 
where 



B = 




(10) 



D is the range interval and E, which will be termed 
the effective surface impedance, is given by 



E{ND,ID) = [ $ + — j exp(-ft) 
V 



(11) 



where ~q has the value appropriate to range {ID) and ^ and 
| (see Fig. 1 and Equation (5)) depend on both ND and ID. 
(In the special case of a plane Earth the right-hand side of 



Equation 
ID only.) 
the factor 



11) simplifies to 7?/r? ; E is then a function of 
C{N,I) is a coefficient which takes into account 



R 



in the integrand. 



y/rlR-r) 
Equation (9) may be written in the form 
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3- 



iV— 1 



G(ND) = 1 -BD % Yj C{N,I)E{ND,ID)G(ID} - BD V2 CiN,N)EIJVd r ND)G{ND) 

7=0 



It will be seen that the term G{ND) appears on both sides 
of this equation. Re-arrangement gives 



N-1 



G(ND)=- 



1 -BD % Yj C{N,I)E{ND,ID)G(ID) 
/= 

1 +BD % C(N,N)E(NDJJD) 



(12) 



G(Q) is then put equal to unity, and successive values of G 
are obtained by repeated applications of Equation {12). 

The problem of choosing the coefficients C is com- 
plicated by the fact that they have to take account of the 
singularity of the kernel of the integral (in Equation (6)) 
at both ends of the range of integration; it is also necessary 
to consider the behaviour of G(r) at the origin. Provided 
that there are no marked irregularities at extremely short 
ranges, it is clear that near the origin G(r) must correspond 
to the groundwave attenuation factor for a uniform plane 
Earth, and it is well known that this may be expressed as a 
power series in r v * . 

The problem is to evaluate 



JE(R,r)G(r) .dr as D Y2 Y" C(N,I)E(ND,ID)G(ID) 

o 

Let f(m) be a function of an integer m such that 

/(0) = 1 | 

f(m) =m- V2 (m>Q) (13) 

It is helpful to factor out this function and write 



The coefficients K(N,I) are chosen on the following 



basis 



N= 1,2or3 

K(NJ) is chosen to perform the integration exactly 
if the integrand in Equation (6), apart from the factor 
jr(R - r)y~ 2 , is a power series in r A . (Note that Sommer- 
feld's groundwave attenuation factor, though singular at the 
origin, can be represented by a power series in r /z .) In 
other words it is assumed that 



N 
E(R,r}G(r) = V* 



<V* 



12 



(16) 



The 7V+1 coefficients K{N,I) are chosen to perform the 
integration exactly for any values of the 7V+1 coefficients 
a n . This is the best that can be done on the assumption 
that E(R,r) is either constant or smoothly varying with r. 

7V=4 



The integral from to D is evaluated from the first three 
ordinates on the assumption that over the first two intervals 



E(R,r)G(r) 



VF 



z> +v ,/2 +V 



(17) 



The integral from 3D to AD is evaluated from the last three 
ordinates on the assumption that over the last two intervals 



E(R,r)G(r) 



■■c + Cl (R-r}+c 2 {R-~r) 7 



(18) 



f(I)f(N-I) 
C(NJ) = — . K(N J) 



fW 



(14) 



so that the sum becomes 



% N 



rrr H fVWN-i)KW)GUD)Ew>jD) 

/uv) / = o 
Equation (12) becomes 



JV-1 
1 -BN V2 D % J2 f(I)f(N~I)K(N,I)G(ID)EWD,ID) 

gind) = Lz ~: UW 



1 +BD Y2 KiN,N}E(NDJVD) 



The factors f(N~I), f(I) in the sum represent the factors 
r~ v \ (R - r)~ " 
last ordinates. 



r , (R - r) 2 in the integral, except for the first and 



The integral from D to 3D is evaluated from the second, 
third and fourth ordinates on the assumption that the 
whole of the integrand is parabolic. 

7V>4 



The integral is evaluated over the first and last intervals as 
for TV = 4. Over the remaining intervals it is evaluated from 
all ordinates except the first and last by means of the 
simplest version of Gregory's formula, using only first dif- 
ferences. This version of Gregory's formula has weighting 

coefficients equal to 5 /l2, 1 ^/12, 1, , 1, ' S /12, 

/12. (It is advantageous to have unit weights for ail 
except a few ordinates at each end, particularly if there are 
discontinuities in the slope or properties of the ground.) 

Table 1 contains the coefficients K{NJ). 

For 7V>5 the coefficients are the same as for N=b with 
N—b unit coefficients inserted between the third and fourth 
coefficients. 
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TABLE 1 

The coefficients K{NJ) 



\ I 














\ 





1 


2 


3 


4 


5 


N \ 














1 


M4159 


2-00000 










2 


0-75605 


2-07238 


0-75605 








3 


0-65170 


2-73801 


-0-40802 


1 -93606 






4 


0-76430 


2-13807 


0-34003 


1 -26667 


1 -20000 




5 


0-76430 


2-22140 


0-27860 


0-89477 


1 -35000 


1-20000 



4. Organisation of program 

4.1. General 

Section 8 is a brief specification of the program, with 
specimens of the data and output. 

The complex groundwave attenuation factor is com- 
puted at a series of ranges nD, where < n <7V. It accepts 
as data heights and ground constants at any number of 
ranges, which need not be regularly spaced, along a radial. 

Heights are specified relative to a spherical datum, 
the radius of which is chosen by the user to take account 
of refraction; alternatively departures from a smooth 
sphere may be ignored. The choice of the factor by 
which the Earth's radius should be multiplied has been 
studied by Rotherham. 1 1 His results suggest factors of 
1-18 and 1-25 for the l.f. and m.f. bands respectively. 
Provision is also made for calculations relating to a smooth 
and plane but inhomogeneous Earth; these are more rapid. 

4.2. Plane-Earth calculation 

It is convenient to begin by considering this simple 
case, which is based on an integral equation due to Wait, 
of which Equation (6) is a generalisation. 

The ground constants are read, in terms either of 
dielectric constant and conductivity, (k and a) or of com- 
plex relative surface impedance (r?/7? ). By interpolation 
and extrapolation an array of values is formed for ranges 
equal to multiples of the calculation interval D. Equa- 
tions (7) and (8) are then used to form a second array, so 
that t?/t? (complex), k and a are stored for each range. 
(Note that in general the results will not be quite the same 
for the two methods of specifying ground constants, since 
conversion from one to the other takes place after inter- 
polation.) 

At the beginning of each run a table of values of/ 
according to Equation (13) is stored and used throughout 
the run without re-calculation. The coefficients K{N,I), 
apart from those that are unity, are stored with the program. 



G(Q) is then set equal to unity and Equation (15) is 
applied repeatedly to store successive values of G(ND). If, 
as is usually the case, N is large, the total number of 
additions required to form the sum N times is approxi- 
mately equal to N 12. Bearing in mind that for a plane 
Earth E{ND,ID) is simply read from the stored array of 
(complex) values of n/v , and that most of the coefficients 
K are unity, it is found that about 1N 2 I2 real multiplica- 
tions are required. 

4.3. Calculation for spherical or irregular Earth 

When the Earth is not treated as plane, the quantity 
E(ND,ID) obtained from Equations (11) and (5) must be 
computed afresh every time the equation is applied to 
compute an additional value of G. 

If heights are read, so that the Earth is not to be 
treated as a smooth sphere, an array of these is formed by 

interpolation and extrapolation for ranges — D, 0,D, , 

ND, and is stored. Next, heights above the spherical 
datum are converted to heights above the tangent plane at 
the transmitter, which is assumed to be parallel to a line 
meeting the surface at ranges ±D. The rest of the calcula- 
tion is the same as that for a plane Earth, except that for 
each step it is necessary to calculate an array of values of E 
from Equation (11). For example, if G(nD) has already 
been calculated for 1 < n <N, and GUVD) is about to be 
calculated, an array of values of E(ND,ID) (0</<JV) is 
first formed and stored. This process occupies almost half 
the running time, even though a rapid method is used for 
computing exp (—/£•) for a series of slowly-changing values 
of?. 

A number of approximations are made: 

(i) No distinction is made between an angle to the hori- 
zontal and its sine or tangent. 

(ii) The spherical datum is replaced by a paraboloid of 
the same curvature at the origin. 

(iii) Distances are regarded as the same whether they are 
measured along a straight line at any inclination or 
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along the spherical datum, except in the computation 
of|. 

iv) In computing % the approximation used is equivalent 



g § the approxirr 
(1 +x 2 / 2 by 1 + 



to replacing (1 +x ) by 1 +x I whenx is small. 



5. Results 

5.1. Uniform plane Earth 

A number of calculations were performed for a uni- 
form plane Earth, with various values of the ground con- 
stants, because values of the complex attenuation factor 
had been published by Norton (see References 4 and 6 for 
the magnitude and phase respectively). 

Norton used as independent variable a complex 
quantity 



Pi = pe -i b ^-y^ Q [— ) r 



(19) 



and termed p the 'numerical distance'.* 

Table 2 shows some results computed for values of 
the ground constants chosen to facilitate comparison with 
Norton's results for b = 30°. For the shortest calculation 
interval (the quantity!) in Equation (13)) used, 250 m, the 
discrepancy with Norton's results is less than 1% in magni- 
tude and (usually) less than 1° in phase. It is believed that 
a few discrepancies in phase greater than 1 can be attri- 
buted to small imperfections in the plotting of Norton's 
curves. An additional means of checking the phase is pro- 
vided by the fact that the attenuation factor should tend 
asymptotically to 180° —b, in this case 150°. Examination 
of the last few entries in the fourth column of the table 
shows that the departure of the phase from 150° decreases 
as the inverse square of the range, with an accuracy of 
about 0-01°. 

The right-hand part of Table 2 shows the effect of 
increasing the calculation interval. The values relative to 
those for D = 250 m may be taken as indicating the errors 
resulting from numerical integration with a finite calculation 
interval when there are no discontinuities in the ground 
constants. It will be seen that these errors tend to be 
greatest for the first value computed, when integration has 
taken place over only one interval, but that after a few 
intervals the errors tend to settle down. 

Referring back to Equation (6), it will be recalled that 
the integral is proportional to the departure from unity of 
the attenuation factor. It follows that when the attenua- 
tion factor is small it is derived as the difference between 
two quantities that are nearly equal. Bearing this in mind 
it seems surprising that the error does not increase markedly 
towards the end of the range. 



* in Norton's papers / is replaced by — : since, unlike the pn 
author, he used exp(-iwt) as time factor. 



5.2. Uniform spherical Earth 

Computations for a uniform spherical Earth test the 
validity of approximations made in the derivations of 
Equation (6), particularly those made on the assumption 
that the gradient of the ground is small. 

Results have been compared with curves published 
by the CCIR, but these curves were calculated manually, 
and can be up to 3 dB in error. More recently the CCIR 
Secretariat has written a computer program for more 
precise calculations, using formulae published by Bremmer; 
curves for one case (dielectric constant 4, conductivity 
3-10~ 3 S/m) have been published. 1 In this case, using a 
calculation interval (£>) of 2500 m, the agreement was to 
the accuracy with which the curves could be read: 0-5 dB. 
The maximum range was 900 km, at which the gradient 
relative to the tangent plane was 0-14. 

In order to achieve a given accuracy for a given 
magnitude of attenuation factor it was found that the 
calculation interval had to be smaller for a spherical Earth 
than for a plane Earth. The reason appears to be associated 
with the factor exp(— j%) in Equation (5), which causes the 
phase of the integrand to change rapidly with r at the 
longer ranges. Large errors can occur when | changes by 
more than 20 from one ordinate to the next, it appears 
advisable to avoid a phase change of more than about 10°; 
the condition for this is 



£><■ 



2a 2 X 
9R 2 



(20) 



In practice this overriding restriction only comes into play 
at ranges of at least several hundred kilometres. 

Table 5 in the Appendix indicates the ranges, and 
attenuation factors, at which specified errors were found to 
occur in practice. 

5.3. Two-part Earth 

Some calculations have been made for a plane Earth 
where the conductivity is finite but uniform along part of 
the path of propagation and infinite along the remainder. 
The results were in good agreement with those published by 
Wait. In general the accuracy for a given attenuation 

factor and number of ordinates was much the same as for a 
uniform plane Earth. Greater errors tended to occur in 
the immediate neighbourhood of the discontinuity, but 
these died out within a few calculation intervals. This 
effect is also apparent in the more complex case to be con- 
sidered in the next section. 

5.4. Test based on the Brewster angle 

In order to test the method of calculation further it 
has been applied to a case in which the ground is neither 
homogeneous nor flat, but in which it is possible to deter- 
mine the attenuation factor by other means. The case 
used for this purpose is designed to exploit the fact that 
vertically-polarized plane waves incident on a perfect 
dielectric are transmitted without reflection when the tan- 
gent of the angle of incidence is equal to the refractive 
index. This angle of incidence is termed the Brewster 
angle. 
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perfectly-conducting 
Case A 



surface impedance G"05 




perfectly-conducting 
Case B 



surface impedance 0-05 




surface impedance 0-05 perfectly-conducting 
Case C 

Fig. 2 - Test based on Brewster angle — sections 

Fig. 2 illustrates three cases for which the ground- 
wave attenuation factor is plotted in Fig. 3(a). The calcu- 
lation interval is 0-5 km, so that the integral is calculated 
from 201 ordinates for a range of 100 km. The frequency 
is 1 MHz. 

In Case A the ground is plane and is perfectly con- 
ducting for the first 50 km from the transmitter. There- 
after it is a non-conducting dielectric having a refractive 
index of 20, corresponding to a (real) relative surface 
impedance equal to 0-05. Fig. 3(a) shows that the attenu- 
ation factor is as expected unity for the first 50 km, de- 
creasing sharply thereafter. 

Case B differs from Case A in that the ground beyond 
50 km curves upwards in such a manner that a radial from 
the transmitter would meet it at the Brewster angle. Since 
there should be no reflection at the surface of the dielectric, 
spherical waves radiated from the transmitter should pass 
across the boundary undisturbed just as if the ground were 
plane and perfectly conducting over the whole of the path 
of propagation. Strictly, the Brewster phenomenon applies 
only to plane waves incident on a plane surface but, in 



view of the smaii curvatures involved, it seems reasonable 
to assume that it will be applicable in the case under con- 
sideration. 

Fig. 3(a) shows that, as expected, the attenuation 
factor remains unity apart from one point at the boundary. 
(It was not obvious what value should be assigned to the 
surface impedance at the boundary; a value of 0-025 was 
chosen.) Fig. 3(b), which is the same as Fig. 3(a) with the 
vertical scale expanded ten times, shows that beyond the 
boundary the attenuation factor is within 0-15% of unity. 

Case C is the same as Case B but with the transmitter 
placed on the curved dielectric surface 50 km from the 
boundary. In this case the attenuation factor is not known 
a priori except for a range of 100 km; at this range recipro- 
city demands that the attenuation factor must be unity to 
agree with Case B. Fig. 3(a) shows that the rapid attenu- 
ation over the dielectric surface is followed by a recovery as 
expected, and that the calculated attenuation factor at 100 
km is unity. (The tabulated results gave it as unity to 
within 0-01%.) 

The results described above provide a useful test of 
the approximations made in dealing with a curved surface. 
in this connection it is worth noting that the maximum 
gradient of the surface in Case B is 0-084. 

The results shown in Figs. 3(a) and 3(b) and discussed 
above were computed for a calculation interval of 0-5 km. 
Fig. 3(c) shows the change in the computed attenuation 
factor for Case C when the calculation interval is increased 
to 4 km (reducing the computation time on a UNIVAC 
1 108 computer from 4 s to 0-1 s for 100 km). The change 
plotted may be regarded as an indication of the error 
resulting from the method of integration used, with ordi- 
nates spaced at 4 km. It will be noticed that this error is 
quite small except for the inevitable disturbance near the 
boundary. (Note that in this case there was no ordinate at 
the boundary.) 

5.5. Calculations over a real path 

Unfortunately accurate information on the ground 
constants over a real path is lacking, so that a comparison 
between theory and experiment was not possible. However 
if it is assumed that tests such as those described above give 
sufficient confidence in the method of calculation used, 
calculations over a real path provide information on the 
importance of the ground contours. 

Fig. 4(a) shows the ground section (neglecting the 
Earth's curvature for convenience in drawing the figure) 
along a radial 300 km long extending from a transmitter at 
Crowborough in Sussex towards Paris. It will be seen that 
to simplify the work of data preparation this profile has 
been greatly simplified. Although in fact the path en- 
counters a cliff at the French coast, the program can only 
take the height to have changed by the appropriate amount 
(150 m) from one ordinate to the next. The cliff is there- 
fore treated as if it were a slope that is inversely propor- 
tional to the calculation interval. For the smaller of the 
two intervals considered, 2 km, the maximum gradient 
supposed by the program to exist is 0-08. 
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case B 




50 
range, km 

Fig. 3 - Test based on Brewster angle — results 

(Frequency 1 MHz) 

(a) Linear plot of attenuation factor (b) As (a) with expanded vertical scale 

(c) Change in calculated attenuation factor resulting from increase in calculation interval from 0-5 km to 4 km 
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effect of increasing D 
from 2 km to 4km 



(d) 



100 



range, km 200 



300 



Fig. 4 - Radial from Crowhorough towards Paris 
(Frequency 1 MHz) 



(a) Ground profile (for convenience shown relative to flat Earth, although Earth curvature was taken into account) 
(b) Attenuation factor (decibels) for a calculation interval of 2 km (taking account of ground profile) 
(c) Effect on attenuation factor of departures from a smooth sphere 
(d) Effect on calculated attenuation factor (taking account of ground profile) of increasing the calculation interval from 2 km to 4 km 
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It must be recognised that the simplifying assumption 
made would be invalidated by the presence of a cliff. It is 
not known what the effect of this would be in practice. 

Atmospheric refraction was assumed effectively to 
increase the radius of the Earth in the ratio 1-25:1. 

The ground on each side of the Channel was taken to 
have a dielectric constant of 4 and a conductivity of 
8 x 10~~ 3 S/m. The sea was assumed to have a dielectric 
constant of 81 and a conductivity of 4-6 S/m. The fre- 
quency was 0-95 MHz (intermediate between the fre- 
quencies of two transmissions that were of interest). 

Fig. 4(b) shows the computed attenuation factor, 
and (c) shows the effect on this attenuation factor of the 
departure of the surface of the Earth from a smooth 
sphere, which was determined by repeating the run with all 
heights set to zero. Fig. 4(d) shows the effect of increasing 
the calculation interval from 2 km to 4 km. 

St will be seen that the effect of the ground contours 
is not great, except for a local increase of about 3 dB in the 
attenuation factor at the top of the steep slope on the 
French side of the Channel. Even in the vicinity of the 
sudden change in height and conductivity the effect of 
doubling the calculation interval is less than 1 dB; it is 
therefore believed that the sharp rise shown in Fig. 4(c) is 
real and not merely an error resulting from the use of a 
finite calculation interval. 

It may be concluded that if mountainous areas are 
excluded the effect of the ground contours is not usually 
important. Apart from one point on Fig. 4(c) the fluc- 
tuations caused by ground contours are less than those 
resulting from the presence of buried conductors and other 
irregularities encountered when making field-strength 
measurements. Nevertheless, the siting of the Crow- 
borough transmitter on high ground gives a general increase 
in the groundwave field strength of the order of 0-8 dB. It 
is true that the subjective effect of such a small improve- 
ment is hardly perceptible to the listener, but the increased 
capital and running costs of a 600 kW transmitter compared 
with a 500 kW transmitter are by no means negligible. 

An even greater increase in groundwave field strength, 
equivalent to a doubling of transmitter power, could have 
been achieved by siting the transmitter at the top of a 
steep slope 150 m high (for example if a transmitter were 
to be sited on the French coast to serve England). 

Reference to Fig. 4(d) shows that a calculation inter- 
val of 4 km is not too great. Increasing the interval from 
2 km to 4 km reduces the calculation time from 2-3 s to 
0-6 s. 

6. Conclusions 

From such tests as have been performed it may be con- 
cluded that the method of computation adopted is free 
from significant errors, provided that the calculation inter- 
val is made sufficiently small. An indication of the maxi- 
mum interval may be obtained from the examples given in 
Section 5; it is discussed in more detail in Section 8. 



Although the computer program was intended pri- 
marily for research purposes, it does appear that the 
routines used are fast enough to be incorporated in a 
program for plotting the ground contours of proposed 
transmitting stations. The principal obstacle to be over- 
come in the preparation of such a program is that of 
obtaining sufficiently detailed information on the ground 
constants. It is difficult to deduce this information from 
field-strength measurements without making assumptions, 
on inadequate evidence, about the variation of the dielectric 
constant with the conductivity. On the other hand it 
would be possible, though not easy, to make measurements 
of phase as well as amplitude and so determine the conduc- 
tivity and dielectric constant independently. 

The computation routines are fast enough for pre- 
dicting service areas, and they provide quite an economical 
way of calculating propagation curves over a uniform 
smooth Earth. Nevertheless the running costs might be- 
come rather high if very large numbers of long paths had to 
be calculated in order to plot contours taking into account 
groundwave interference. 

If the rather crude integration procedures were to be 
refined it is possible that calculation intervals could be 
doubled, reducing running times by a factor of 4. 
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Appendix 
Specification of Program GW1 



8.1. Outline of program 

The groundwave attenuation factor is tabulated as a 
function of range taking into account an arbitrary variation 
of the ground constants along the path. The Earth may be 
plane, in which case the program is faster, or have an 
arbitrary variation of height, relative to a spherical datum, 
along the path. The table may show a comparison of the 
attenuation factor with that computed for another case 
dealt with earlier in the same run. The program is written 
in FORTRAN 5 for a UNIVAC 1108 computer. 

The calculation is performed by two subroutines, 
INT1 (plane Earth) and INT3 (curved Earth). 



GW1 (main 
program calls: 



INT1 calls 
INT3 calls 



PRODATto take in height and ground- 
constant data. An additional 
entry PROF, called by GW1, 
forms it into arrays by linear 
interpolation. 

INT1 to form an array containing 

the (complex) attenuation fac- 
tor for a plane Earth at equal 
intervals of range. There are 
additional entry points INT1A, 
called by GW1, and INT, not 
called by GW1. 

INT3 As INT1, but for a curved 
Earth. There is an additional 
entry point INT3A called by 
GW1 

POL for converting complex num- 

bers to polar form, with (op- 
tional) magnitude also in DB. 
There is provision for zero 
argument, 
a special integration routine 



INT2 

INT2 

CIS 



a time-saving exp(pc) function. 



For each case the attenuation factor is calculated over the 
range specified and stored in an array before the tabulation 
of results begins. 

8.2. Facilities available 

Terrain data in the form of heights (optional) and 
ground constants may be specified for one range or for a 
series of ranges, which need not be in arithmetical pro- 
gression, and need not correspond to ranges at which the 
attenuation factor is to be computed. At intermediate 
ranges values are computed by linear interpolation. Be- 
yond the extreme ranges for which data are supplied, the 
height, dielectric constant and conductivity will be regarded 
as constant. For example, if terrain data are provided at 
three ranges, and if the heights at these ranges are repre- 
sented by the points plotted in Fig. 5, height is assumed to 



range 

Fig. 5 - Interpolation and extrapolation of terrain data 

vary with range according to the continuous line. Data for 
one range are sufficient to specify constant height and 
ground constants. 

If no heights are specified the Earth is assumed to be 
plane or spherical. 

Ground constants may be specified in terms of relative 
surface impedance, assumed to be equal to (e r + 1) 
where e r is the complex relative permittivity, or in terms of 
real relative permittivity and conductivity. Whichever 
form is specified, both are tabulated. 

In general, heights are specified relative to a spherical 
datum of radius 6371i? F km where i? F is a factor taking 
account of refraction. Rotherham's results suggest that 
1-25 should be used for the m.f. band and 1-18 for the l.f. 
band. 

I NTS may be used for a plane Earth, or for heights 
specified relative to a plane datum, by specifying a large 
value of i? F (e.g. 1000000). Alternatively INT1, which is 
faster, may be used for a plane Earth. 

The same terrain data may be used for two or more 
successive cases. For example, it may be desired to re- 
compute the results for a new calculation interval, fre- 
quency, or datum radius. Having provided both height 
and ground-constant data, it is possible to ignore the 
heights and re-calculate for a plane Earth. (But to re- 
calculate for a smooth spherical Earth it is necessary to re- 
supply the ground-constant data.) 

The tabulation interval (in terms of range) may be 
chosen to be equal to the calculation interval or to any 
integral multiple of it. 

In addition to tabulating the attenuation factor for 
each case computed, the table may include a comparison 
with a case dealt with earlier in the same run. This facility 
is provided by using two alternative sets of locations in 
which the attenuation factor for each case may be stored. 
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If a comparison is called for the results are compared with 
the only previous results available, i.e. those for the last 
case stored in the other set of locations. 

Note that the attenuation factor is stored for all the 
values calculated (at the calculation interval) not merely 
for those tabulated. Thus if for Case 1 the attenuation 
factor is calculated at 1 km intervals and tabulated at 5 km 
intervals, and if for Case 2 it is calculated at 2 km intervals 
and tabulated at 4 km intervals, the results for Case 2 may 
include a comparison with Case 1 in every line of the table. 

Nevertheless a comparison may be called for when 
values of the attenuation factor for the previous case are 
not all available. For example, if Case 1 calls for calcula- 
tion of ranges of (2) 100 km, and Case 2 at (5) 200 km 
(with tabulation at all these ranges), and if a comparison is 
called for, comparison will be possible only at ranges (10) 
100 km. This is permissible; the comparison figures will 
be given only on those lines of the Case 2 table correspond- 
ing to ranges for which results were stored in Case 1. Two 
ranges are regarded as equal when they differ by less than 
one part in 10 5 . 

8.3. Data 

All data are numerical, and in free-field format, each 
number being terminated by a comma. Floating-point 
representation is not permitted. Symbols in Roman type 
below are FORTRAN variables. Their units differ from 
those used in the main text of this report. 

For each case the following cards are required: 

First Card 

NCASE, 

This is an integer for which (KNCASE<10QG0. It 
represents the case number, and is printed with the results. 

Second Card 

RF, FQ, D, Nl, LOG, NT, NCOMP, M, 

where: 

RF represents the ratio in which the radius of 

the Earth is effectively increased by refrac- 
tion 

If RF is made zero, this is taken as an indi- 
cation that the plane-Earth subroutine INT1 
is to be used. Any heights read with ground- 
constant data will then be disregarded (but 
will nevertheless be stored, so that a subse- 
quent run with the same terrain data and a 
non-zero value of RF will take account of 
them) 

FQ is the frequency in MHz 

D is the calculation interval (also the inte- 

gration interval) in km 

Nl an integer, is the maximum range as a 

multiple of D. It should preferably be an 



integral multiple of NT {see below). Nl 
must not exceed NIM, a parameter of GW1 
at present set at 500 

LOG an integer, must be either 1 or 2. It indi- 
cates which set of locations is to be used for 
storage, with a view to comparison between 
cases (see NCOMP below) 

NT is the tabulation interval as a function of D. 

Results will be printed for ranges 0, NT.D, 

2NT.D, , m.NT.D where m.NT is the 

largest multiple of NT not exceeding Nl 

NCOMP an integer, must be either zero, if no com- 
parison is required, or 1 if a comparison is 
required. In the latter case the current 
values of the attenuation factor will be com- 
pared with those stored in the other set of 
locations 

M is an integer such that 0<M<4. It has the 

following significance 

Previous terrain data to be used. The 
next card to be read will be the next 
value of NCASE. 

1 Heights all zero. Permittivity and con- 
ductivity to be read. 

2 Height, permittivity and conductivity 
to be read. 

3 Heights all zero. Relative surface 
impedance to be read. 

4 Heights and relative surface impedance 
to be read. 

Third and subsequent cards (Terrain data) 

These, like other data, are in free-field format, each number 
being terminated by a comma. The way in which the 
information is interpreted is outlined in Section 8.2 (above! 
and Fig. 5. There may not be more than NCM cards (in- 
cluding the terminating card). NCM is a parameter of GW1, 
at present set to 200. 



Each card carries 3 numbers (M 
numbers (M = 2 or M = 4). 



1 or M = 3) or 4 



M = 0: No terrain data will be supplied. The terrain data 
last provided will have been stored as it was read 
from cards, and this will be re-processed for the 
new frequency, calculation interval, etc. 

M = 1 : Range,k,a, 

The range is in km. The range must increase 
from card to card, apart from the terminating card 
k is the dielectric constant, i.e. the real part of the 
relative permittivity 
a is the conductivity in Siemens/m 

M = 2: Range,height,£,a, 

Heights, in feet, are relative to the spherical datum 
of radius specified. The curved-Earth subroutine 
requires a height at a range — D, where D is the 
calculation interval, since heights at — D and +D 
are used by the subroutines to derive the gradient 
at range zero. It is important to bear this in mind 
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when preparing terrain data, and to provide at 
least one height at a range <— D if the ground 
slopes near the transmitter.* 

M = 3: Range,real part of ETAR,imaginary part of ETAR, 
ETAR = t?/t? is the relative surface impedance, 
which is taken to be equal to\k + 1 -/(18000a/ 
FQ)f- y ' 

M = 4: Range,height,real part of ETAR, imaginary part of 
ETAR, 



TABLE 3 
Specimen set of data for six cases 



Terrain data terminating card 

For any value of M>0 the terrain data must be 
followed by a terminating card (which will be followed by 
the next case number). This should contain at least as 
many numbers as each of the terrain data cards. The first 
number must be less than the range on the preceding card. 
This reduction in the range is taken to indicate that the 
terrain data is complete. The other numbers on the card 
will be ignored, but they must be there (otherwise the 
program will read other cards to find them and the remain- 
der of the run will produce on|y rubbish*). The terminat- 
ing card may contain 0,0,0,0, 

Table 3 is a specimen set of data for six cases, 
numbered 21 to 26. All are for a spherical datum (RF = 
1-25), a frequency of 0-95 MHz, and a range of 300 km 
(304 km for Case 24). The path includes some sea. 

In Case 21, heights, relative permittivity and con- 
ductivity are given at 15 points. The attenuation factor 
is to be calculated at intervals of 2 km and tabulated every 
4 km. The results are to be stored in Locations 1. 



Cases 22 to 25 are the same, except that the attenua- 
tion factor is to be calculated and tabulated every 4, 6, 8 
and 10 km respectively. The results are to be stored in 
Locations 2 and compared with Case 21 (in order to see 
how the errors increase as the calculation interval is 
increased). 

Case 26 is the same as Case 21 except that all heights 
are zero. The results are to be stored in Locations 2 and 
compared with Case 21 (to see how much effect the ground 
contours have). 

8.4. Output 

Table 4 is the output for Case 25, the data for which 
is included in Table 3. 

Separate running times are printed for 'data', 'calcu- 
lation' and 'table'. 'Data' includes input together with 
processing of input to produce arrays containing heights 



' Dr. J.H. Causebrook has recently pointed out that a height for a 
negative value of range is not strictly necessary, since the assign- 
ment of the wrong gradient to the tangent plane has no effect on 
the result. 

'* If and when a new version of the program is produced, the input 
arrangements should be changed to avoid this pitfall. 



(above a spherical datum) and relative surface impedance, 
both at multiples of the calculation interval. 

'Calculation' ends with the attenuation factor stored 
in a complex array for a set of ranges equal to multiples of 
the calculation interval. 

'Tabulation' includes ancillary calculations such as 
conversion of complex numbers to polar form and compari- 
son between cases. 

The bulk of the table is self-explanatory except for 
the following: 

ATTEN(G) is the complex groundwave attenuation 
factor, i.e. the complex ratio of the field strength to that 
which would exist if the Earth were plane and perfectly 
conducting with a path length equal to the straight line 
(chord) between transmitter and receiver. 

Phases are in degrees. 

In the comparison columns (omitted if NCOMP is 

zero) 

'CASES (25 - 21)' means 'G(Case 25) - G(Case 21)' 

'CASES (25/21)' means 'G(Case 25)/G(Case 21)'. 



(RA-172) 



15 



I f i 1 1 



• I • 



i • I I 



I J i I i I ! 1 1 ( » » i 1 



^ eg 

5 S 

a- 



< ■ - < ; • * ; •- 



- i I ! ! ! I I ( I 1 I i i >, t i I I i < - ( I I '. ( < i 5 ( 



•: : i < i 



(RA-172) 



16 



The first line of the table (omitted in plane-Earth 
cases) contains only a negative range and a height. The 
range is — D, where D is the calculation interval (not the 
tabulation interval if this is different) and the height is that 
interpolated at that range. See Section 8.3, under Terrain 
data (M = 2)' for the significance of this. 

8.5. Error indications 

Sufficient data is reproduced in the output to enable 
errors in data to be diagnosed. However there are no 
checks for illegal data, and the use of free-field format 
(which is regretted) means that the omission of numbers 
can produce rubbish in the remainder of the run. Very 
long runs are therefore inadvisable. 

8.6. Accuracy 

Known sources of error are as follows. 

8.6.1. Integral equations 

Both the integral equation used by INT1 and the 
more general integral equation used by INT3 are based on 
an approximate method of reducing a two-dimensional 
equation to a one-dimensional equation. Such checks as 
have been made (see Reference 7, pp. 109— 110) indicate 
that the error does not exceed 0-1 dB in practice. 

8.6.2. Finite gradients 

The curved-Earth routine I NTS contains approxi- 
mations based on the assumption that the gradient of the 
ground, relative to a plane datum, is small. The only 
evidence as to the effect of this is as follows. 

In the case of a prefectly-conducting spherical Earth 
at 1 MHz, calculation out to a range of 900 km gives agree- 
ment with published curves (which cannot be relied upon 
to less than 0-5 dB). The maximum gradient relative to 
the tangent plane was about 0-14 in this case. 

An artificial case, based on the fact that waves incident 
on a pure dielectric at the Brewster angle are wholly 
absorbed, involved a gradient increasing from 0-05 to 0-085 
over the path. The error did not exceed 0-01 dB. 



However, in general there may be an error equal to the 
effect of moving the discontinuity by a distance up to half 
the calculation interval. For example, in Case 25 (see 
Table 4 for output), which uses the same terrain data as 
Case 21, the calculation interval is 10 km. Identical 
results would have been obtained with the boundary any- 
where between 30 and 40 km. 

The last column of Table 4 indicates that, with a SOX 
calculation interval, errors between 1 and 2 dB can occur 
in the neighbourhood of discontinuities. These errors 
diminish roughly in inverse proportion to the interval, and 
they tend to die away at points remote from the discon- 
tinuities causing them, unless a discontinuity is near to the 
transmitter. 

If accurate results are required in the immedaite 
neighbourhood of discontinuities, or if there is a discon- 
tinuity near to the transmitter, a small calculation interval 
(e.g. 2 km in the m.f. band) cannot be avoided. In other 
cases it is recommended that the data should smooth out 
discontinuities. This may be done by providing data cards 
at ranges R 1 ± VtD where R t is the range of the discon- 
tinuity and D the calculation interval. 

8.6.4. Effect of finite calculation interval in the 
absence of discontinuities 

Even in the absence of discontinuities error is 
caused by the finite calculation interval. If ranges less than 
about ED are discounted, it is found that the error (in dB) 
increases as the magnitude of the attenuation factor 
decreases. 

Plane Earth 



Sufficient results have not been obtained to specify 
the maximum error with confidence, but the following 
estimates are based on such results as are available.* 

For an error less than 0-1 dB, the minimum number 
of calculation intervals necessary is as follows: 

Attenuation factor at 

extreme range (dB) -20 -30 -40 

Number of intervals 30 50 170 



It may be concluded that the error due to finite 
gradients will rarely be important in practice. 

8.6.3. Neighbourhood of discontinuities 

Since terrain data is linearly interpolated discon- 
tinuities are not, strictly speaking, possible. However the 
effect of a discontinuity may be obtained — e.g. at a land/ 
sea boundary — by providing different data for two ranges 
that are close together. See, for example. Table 3, Case 21, 
ranges 3T5 and 32-5 km, representing a land/sea boundary 
at 32 km. Now the program takes account of the height 
and ground constants only at ranges equal to multiples of 
the calculation interval. The most accurate result will be 
obtained when the discontinuity happens to be half-way 
between two ranges for which calculation is performed. 



If larger errors can be accepted the interval can be 
increased roughly in proportion to the acceptable 
error. 

Spherical Earth 



For a given accuracy smaller calculation intervals are 
required for spherical Earth than for a plane Earth. The 



* Since this was written Dr. J.H. Causebrook has reported that the 
minimum numbers of intervals indicated proved inadequate for 
poorly-conducting ground at the higher medium frequencies. It 
is recommended that in case of doubt the calculation should be 
performed twice for different intervals and that the comparison 
facility should be used. 
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TABLE 5 

Errors resulting from finite calculation interval 
(The Earth radius factor was unity in all cases) 



Frequency 
(MHz) 


Dielectric 

Constant 


Conductivity 
(Sim) 


Interval 
(km) 


No. of 
Intervals 


Range 
(km) 


Attenuation Factor 


Error in 
Magnitude (dB) 


Magnitude (dB} 


Phase 


0-2 


4 


0-003 


5 


156 
162 


780 
810 


-28 
-29 


-438° 
-459° 


0-1 
0-2 


10 


62 
78 
81 
85 


620 
780 
806 
850 


-22 
-27 

-29 
-30 


-327° 

-440° 
-456° 

-492° 


0-1 
0-2 

0-5 
1-0 


20 


15 

32 
41 


300 
640 
820 


-10 
-23 
-29 


-175° 
-338° 
-468° 


0-2 
0-5 
1-0 


0-7 


4 


0-01 


3-333 


105 
114 
147 


350 
380 
490 


-31 
-34 
-41 


-317° 
-342° 
-450° 


0-1 
0-2 
0-5 


1-0 


4 


0-003 


1 


144 


144 


-38 


-213° 


0-1 



maximum permissible interval decreases as the magnitude 
of the attenuation factor decreases and as its phase retar- 
dation increases. Moreover significant errors can occur if 
the phase f changes by more than 20° from one ordinate 
to the next. It therefore appears desirable to avoid a phase 
change of more than 10°; the condition for this is 



D<- 



2fl 2 X 
9R 2 



(20) 



Sufficient results have not yet been obtained to give a 
formula relating the error to the .variables affecting it. 
However Table 5 summarises the information obtained by 
comparing results for different calculation intervals. The 
error shown (0-1, 0-2, 0-5 or 1-0 dB) represents the ratio 
of the attenuation factor for the interval stated to that for 
a smaller interval. In each case the range quoted is the 
range (ignoring the first few points calculated) at which 
the stated error is first reached. 



8.7. Running time 
The total running time is made up as follows: 
For curved earth (I NTS): 

seconds 



overhead 

+ 

+ 

+ 

+ 

or 



0-06 

0-002 per terrain card 

0-0008NI 

0-000KNI) 2 

0-005/line of table (without comparison! 

0-010/line of table (with comparison) 



The term in (Nl) 2 is the principal component of the calcu- 
lation time. In the case of a plane-Earth calculation 
(INT1) it reduces to 0-00004 (N I) 2 



SMW/VY 
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